logo of company

Intersections of Adversity & Neurodiversity: Adverse Childhood Experiences’ Association to Mental Health & the Buffering Role of Flourishing

This study extends analyses conducted by Bethell et al. (2023), who explored the role of family resilience and connection in promoting flourishing among U.S. children facing adversity. Our study applies a similar framework but focuses on neurodivergent children, examining how adverse childhood experiences (ACEs) influence mental health outcomes and how flourishing behaviors—such as curiosity, emotional control, and task persistence—may buffer these effects. This study was conducted by the Slopen Laboratory at the Harvard T.H. Chan School of Public Health.

Author: Adrián Medina

Date: October 8, 2024

Project Overview

This repository contains files relevant to a study investigating the relationship between adverse childhood experiences (ACEs) and mental health outcomes in children with neurodevelopmental disorders (NDDs), such as Autism Spectrum Disorder (ASD) and ADHD. The study specifically examines the moderating role of child flourishing behaviors—curiosity in learning, emotional control, and task persistence—on the effects of ACEs on mental health outcomes.

Study Design:
The study utilizes data from the National Survey of Children’s Health (NSCH), focusing on children aged 6-17 with reported neurodevelopmental issues (N = 44,776, MAge = 12.2). This cross-sectional analysis examines ACE exposure and its impact on three primary mental health outcomes: anxiety, depression, and behavioral issues. The moderating role of child flourishing behaviors is assessed using interaction terms in logistic regression models.

Data Collection:
Key variables collected include:

  • Adverse Childhood Experiences (ACEs): Derived from parent-reported experiences such as exposure to violence, household dysfunction, and discrimination.
  • Mental Health Outcomes: Parent-reported diagnoses of anxiety, depression, and behavioral problems, confirmed by health care providers or educators.
  • Child Flourishing Metrics: Based on parent-reported behaviors such as curiosity in learning, emotional control, and task persistence, which were used to develop a Child Flourishing Index (CFI).

Analytic Approach:
The primary analytic approach involved logistic regression models to assess the association between ACEs and mental health outcomes, with interaction terms to explore the moderating effect of child flourishing. Key covariates include age, sex, race/ethnicity, and socioeconomic status.

Goals:
The study aims to:

  • Examine the dose-response relationship between the number of ACEs and risks of mental health challenges in neurodivergent children.
  • Identify the protective role of child flourishing behaviors, offering insights into resilience-promoting interventions for this vulnerable population.
Tip

For detailed information on the dataset variables, refer to the NSCH Data Dictionary.

Code Workflow

  1. Data Frame Initialization
  2. Analytic Data Preparation
  3. Frequencies & Descriptive Statistics
  4. Inferential Statistics
  5. Session Information

Data Frame Initialization

Set up the R environment by configuring the CRAN repository, installing essential packages, setting the base path, and loading the primary dataset into the local environment.

Expand Code 
# Set CRAN repository for consistent package installation
options(repos = c(CRAN = "http://cran.rstudio.com/"))

# Install and load the pacman package for efficient package management
if (!require(pacman)) install.packages("pacman")
library(pacman)

# Use p_load function to install (if necessary) and load packages
p_load(dplyr, tidyr, tidyverse, knitr, ggplot2, Hmisc, broom, stats, kableExtra, webshot2, sjPlot, margins, gtsummary, ggstats, broom.helpers, patchwork, sjlabelled, sjmisc, vtable, DiagrammeR)

# Specify the 'base_path' where you can find your data files, ASSUMING they're in the same directory, & set it as WD
base_path <- "~/GitHub/Adversity-Neurodiversity/data_files"
setwd(base_path)

# Load primary data file
neurodivergent_data <- read_csv("neurodivergent_data.csv")
1
See details under the Data Extraction & Filtering (Archived Code) callout below!
Data Extraction & Filtering (Archived Code)

The file being used for these analyses is a subset of a “master” file, Data2_2016to2022.csv, containing the compiled contents of the data releases from the National Survey on Child Health’s 2016-2022 cycles. The master file is nearly 500 MB (0.5 GB) in size, making it cumbersome to process in real-time. For efficiency and ease of analysis, I am using a smaller subset that contains only the relevant data. This subset was generated using the data extraction and filtering process detailed in the archived code below, in an effort to maintain transparency and ensure reproducibility of the workflow.

Expand Code 
# Define the path to the dataset
data_path <- "~/Downloads/Data2_2016to2022.csv"  
Data2_2016to2022 <- read.csv(data_path)

# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status
neurodivergent_data <- Data2_2016to2022 %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis, k2q35a, k2q35c, k2q38a,
         k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, k2q37a, k2q37c, downsyn, downsyn_desc, 
         k2q31a, k2q31c, k2q61a, cerpals_desc, k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, k7q85_r, k7q84_r) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data, "neurodivergent_data.csv")

Analytic Data Preparation

Recode analytic variables and covariates, including ACE counts, Childhood Flourishing Index (CFI), binary outcomes (anxiety, depression, behavioral issues), education, income, sex, and race/ethnicity variables into categorical factors for analysis.

Expand Code 
# Recode 'ace1' to dichotomous variable based on analytic specifications
neurodivergent_data$ace1_recode <- ifelse(neurodivergent_data$ace1 %in% c(2, 3, 4), 1, 
                                          ifelse(neurodivergent_data$ace1 == 1, 0, NA))

# Recode ACE counts and calculate Childhood Flourishing Index (CFI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Total ACE count for each child excluding missing values
    ace_total = rowSums(select(., ace1_recode, ace3:ace12) == 1, na.rm = TRUE),
    
    # Categorize total ACEs for further analysis
    ace_counts = cut(ace_total, breaks = c(-1, 0, 1, 3, Inf), labels = c('0 ACEs', '1 ACE', '2-3 ACEs', '4+ ACEs'), right = TRUE),
    ace_counts = factor(ace_counts, levels = c("0 ACEs", "1 ACE", "2-3 ACEs", "4+ ACEs")), # Recreate 'ace_counts' factor

    # Calculate the Childhood Flourishing Index (CFI)
    CFI = rowSums(select(., k6q71_r, k7q85_r, k7q84_r) == 1, na.rm = TRUE),
    
    # Recode CFI into dichotomous variable
    CFI_dich = case_when(
      CFI == 3 ~ "Flourishing",
      CFI %in% 0:2 ~ "Not Flourishing",
      TRUE ~ NA_character_
    ),
    CFI_dich = factor(CFI_dich, levels = c("Not Flourishing", "Flourishing")) # Recreate 'CFI_dich' factor
  )

# Recode binary outcomes for anxiety, depression, and behavioral issues into factors
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode Anxiety (k2q33b)
    Anxiety = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    
    # Recode Depression (k2q32b)
    Depression = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    
    # Recode Behavioral Problems (k2q34b)
    Behavioral_Problems = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes"))
  )

# Recode education and income categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode education into categorical factors
    highgrade_tvis_cat = case_when(
      higrade_tvis == 1 ~ "Less Than High School",
      higrade_tvis == 2 ~ "High School/Vocational",
      higrade_tvis == 3 ~ "Some College/Associate Degree",
      higrade_tvis == 4 ~ "College Degree/Higher",
      TRUE ~ NA_character_
    ),
    highgrade_tvis_cat = factor(highgrade_tvis_cat, levels = c("Less Than High School", "High School/Vocational", "Some College/Associate Degree", "College Degree/Higher")),

    # Recode federal poverty level categories into factors
    fpl_cat = case_when(
      fpl %in% c(50:99) | fpl_mean < 100 ~ "<100%",
      fpl %in% c(100:199) | fpl_mean < 200 ~ "100%-199%",
      fpl %in% c(200:299) | fpl_mean < 300 ~ "200%-299%",
      fpl %in% c(300:399) | fpl_mean < 400 ~ "300%-399%",
      fpl %in% c(400:999) | fpl_mean < 999 ~ "≥400%",
      TRUE ~ NA_character_
    ),
    fpl_cat = factor(fpl_cat, levels = c("<100%", "100%-199%", "200%-299%", "300%-399%", "≥400%"))
  )

# Recode sex and race categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode sex into categorical factors
    sc_sex_cat = case_when(
      sc_sex == 1 ~ "Male",
      sc_sex == 2 ~ "Female",
      TRUE ~ NA_character_
    ),
    sc_sex_cat = factor(sc_sex_cat, levels = c("Male", "Female")),

    # Recode race and ethnicity
    race = case_when(
      sc_race_r == 1 ~ "White",
      sc_race_r == 2 ~ "Black",
      sc_race_r == 3 ~ "American Indian or Alaska Native",
      sc_race_r %in% 4:5 ~ "Asian, Native Hawaiian, or Pacific Islander",
      sc_race_r %in% 6:7 ~ "Other or Mixed Race",
      TRUE ~ NA_character_
    ),

    ethnicity = case_when(
      sc_hispanic_r == 1 ~ "Hispanic/Latino",
      sc_hispanic_r == 2 ~ "Non-Hispanic/Latino",
      TRUE ~ NA_character_
    ),

    # Recreate 'sc_race_cat' from race and ethnicity combinations
    sc_race_cat = case_when(
      race %in% c("White") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic White",
      race %in% c("Black") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Black or African American",
      race %in% c("American Indian or Alaska Native") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic American Indian or Alaska Native",
      race %in% c("Asian, Native Hawaiian, or Pacific Islander") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander",
      race %in% c("Other or Mixed Race") & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Other or Mixed Race",
      ethnicity == "Hispanic/Latino" ~ "Hispanic/Latino of Any Race",
      TRUE ~ NA_character_
    ),
    sc_race_cat = factor(sc_race_cat, levels = c("Non-Hispanic White", "Non-Hispanic Black or African American", 
                                                "Non-Hispanic American Indian or Alaska Native", "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander", "Non-Hispanic Other or Mixed Race", "Hispanic/Latino of Any Race"))
  )
1
ace1 is originally coded by NSCH as a frequency measure, but it’s needed as indicative variable of presence like the other ace# variables.
2
Part of this includes the calculating of the Childhood Flourishing Index (CFI) based on Bethell et al. (2019) criteria

Frequencies & Descriptive Statistics

Calculate the prevalence of neurodevelopmental disorder (NDD) diagnoses by defining the relevant variables, computing the number of total respondents for each variable, and determining the percentage of diagnoses. Display the results in a table with a table note explaining the presentation format.

Expand Code 
# Calculate prevalence (%) of NDD diagnoses among subject sample
variables <- c("k2q31a", "k2q35a", "k2q38a", "k2q30a", "k2q36a", "k2q60a", "k2q37a", "downsyn", "k2q61a")
variable_names <- c("ADHD", "Autism/ASD", "Tourette’s Syndrome", "Learning Disability", 
                    "Development Delay", "Intellectual Disability", 
                    "Speech/Other Language Disorder", "Down Syndrome", "Cerebral Palsy")

# Pre-calculate total respondents to avoid repeated computation
total_respondents <- colSums(!is.na(neurodivergent_data[variables]))

# Vectorized prevalence calculation
prevalences <- colSums(neurodivergent_data[variables] == 1, na.rm = TRUE) / total_respondents * 100
prevalence_table <- data.frame(Diagnosis = variable_names, Prevalence = prevalences)

# Display prevalence data in a table
prevalence_table %>% 
  knitr::kable("html", caption = "Prevalence of Neurodivergent Conditions") %>% 
  kableExtra::footnote(general = "Prevalence values are presented as percentages for neurodevelopmental disorder diagnosis.")
Prevalence of Neurodivergent Conditions
Diagnosis Prevalence
k2q31a ADHD 58.9175975
k2q35a Autism/ASD 15.2625563
k2q38a Tourette’s Syndrome 1.6302642
k2q30a Learning Disability 39.6132560
k2q36a Development Delay 31.7615419
k2q60a Intellectual Disability 5.6693867
k2q37a Speech/Other Language Disorder 36.2469766
downsyn Down Syndrome 0.9816886
k2q61a Cerebral Palsy 1.4784403
Note:
Prevalence values are presented as percentages for neurodevelopmental disorder diagnosis.

Summarize key demographic characteristics of the neurodivergent sample, including descriptive statistics for age, educational attainment, sex, socioeconomic status, and race. Display the results in a table and visualize distributions using raincloud and bar plots for each demographic variable.

Expand Code 
# Calculate descriptives of age for neurodivergent participants
neurodiv_age_summary <- neurodivergent_data %>%
  dplyr::summarise(
    score_mean = mean(sc_age_years, na.rm = TRUE),  # Mean
    n = n(),  # Sample size
    sem = sd(sc_age_years, na.rm = TRUE) / sqrt(n()),  # Standard error of the mean
    sd = sd(sc_age_years, na.rm = TRUE),  # Standard deviation
    min_age = min(sc_age_years, na.rm = TRUE),  # Minimum
    q1_age = quantile(sc_age_years, 0.25, na.rm = TRUE),  # 1st Quartile
    median_age = median(sc_age_years, na.rm = TRUE),  # Median
    q3_age = quantile(sc_age_years, 0.75, na.rm = TRUE),  # 3rd Quartile
    max_age = max(sc_age_years, na.rm = TRUE),  # Maximum
    IQR_age = quantile(sc_age_years, 0.75, na.rm = TRUE) - quantile(sc_age_years, 0.25, na.rm = TRUE),  # IQR (Q3 - Q1)
    .groups = 'drop'
  )

# Create a neat table using kable
neurodiv_age_summary %>%
  knitr::kable("html", col.names = c("Mean", "Sample Size", "SEM", 
                                    "SD", "Min", "Q1", "Median", "Q3", "Max", "IQR"), caption = "Descriptive Statistics of Age for Neurodivergent Participants")
Descriptive Statistics of Age for Neurodivergent Participants
Mean Sample Size SEM SD Min Q1 Median Q3 Max IQR
12.17523 44776 0.0159532 3.375748 6 9 12 15 17 6
Expand Code 
# Raincloud plot of ages for neurodivergent participants with color
ggplot(neurodivergent_data, aes(x = 1, y = sc_age_years)) + # Set x to 1, as there's only one group
  PupillometryR::geom_flat_violin(adjust = 1.5, trim = FALSE, alpha = 0.5, color = NA, fill = "deeppink", position = position_nudge(x = 0.1, y = 0)) +
  geom_point(alpha = 0.1, position = position_jitter(width = .05), size = .10, shape = 20, color = "deeppink") +
  geom_boxplot(outlier.shape = NA, alpha = 0, width = 0.1, position = position_dodge(width = 0.3), colour = "black") +
  geom_point(data = neurodiv_age_summary, aes(x = 1, y = score_mean), shape = 18, color = "deeppink", size = 3, position = position_nudge(x = 0.1)) +
  geom_errorbar(data = neurodiv_age_summary, aes(x = 1, y = score_mean, ymin = score_mean - sem, ymax = score_mean + sem), width = 0.05, color = "deeppink", position = position_nudge(x = 0.1)) +
  labs(
    title = "Age Distribution of Neurodivergent Participants", 
    y = "Age (Years)", 
    x = "") + # No x-axis label since there's only one group
  scale_y_continuous(breaks = seq(5, 18, by = 1), limits = c(5, 18)) + # Set y-axis range and increments
  theme_minimal() +
  theme(axis.text.x = element_blank()) # Remove x-axis text

Expand Code 
# Bar plot for educational attainment counts categories
ggplot(neurodivergent_data %>% filter(!is.na(highgrade_tvis_cat)), aes(x = highgrade_tvis_cat)) +
  geom_bar(aes(fill = highgrade_tvis_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  labs(title = "Distribution of Adults' Highest Educational Attainment",
       x = "Educational Attainment",
       y = "Count",
       fill = "Educational Attainment") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Expand Code 
# Bar plot for sex counts distribution
ggplot(neurodivergent_data, aes(x = sc_sex_cat)) +
  geom_bar(aes(fill = sc_sex_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  scale_fill_manual(values = c("Male" = "lightblue", "Female" = "pink")) +
  labs(title = "Distribution of Sex",
       x = "Sex",
       y = "Count") +
  theme_minimal()

Expand Code 
# Bar plot for SES counts categories
ggplot(neurodivergent_data, aes(x = fpl_cat)) +
  geom_bar(aes(fill = fpl_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(title = "Distribution of Socioeconomic Status (FPL)",
       x = "Federal Poverty Level Groups",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

Expand Code 
# Bar plot for race counts categories
ggplot(neurodivergent_data, aes(x = sc_race_cat)) +
  geom_bar(aes(fill = sc_race_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust=-0.3, size=3.5) +
  labs(title = "Distribution by Race",
       x = "Race",
       y = "Count") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Calculate and display the number of subjects classified as ‘Flourishing’ based on a Childhood Flourishing Index (CFI) score of 3. Visualize the distribution of children by CFI score using bar plots for both the raw CFI scores and the dichotomized CFI variable.

Expand Code 
# Count the number of subjects with a CFI score of 3 using a logical condition, as this is considered 'Flourishing.'
num_subjects_cfi_3 <- sum(neurodivergent_data$CFI == 3, na.rm = TRUE)
cat("Number of Subjects Flourishing (CFI Score of 3):", num_subjects_cfi_3, "\n")
Number of Subjects Flourishing (CFI Score of 3): 3265
Expand Code 
# Bar plot of CFI counts categories
ggplot(neurodivergent_data, aes(x = factor(CFI), fill = factor(CFI))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "Childhood Flourishing Index (CFI) Score",
       y = "Number of Children",
       fill = "CFI Score",
       title = "Distribution of Children by CFI Score") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

Expand Code 
# Bar plot of CFI_dich counts categories
ggplot(neurodivergent_data, aes(x = CFI_dich)) +
  geom_bar(aes(fill = CFI_dich), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(title = "Distribution of CFI Dichotomous Variable",
       x = "Child Flourishing Index (Dichotomous)",
       y = "Count",
       fill = "CFI") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

Construct and display a frequency table that shows the number of participants with Anxiety, Depression, and Behavioral Issues based on responses indicating “Yes” for each condition. Present the results in an HTML table format using kable.

Expand Code 
# Constructing the frequency table for Anxiety, Depression, and Behavioral Issues where the response is "Yes"
neurodiv_freq_table <- neurodivergent_data %>%
  dplyr::summarise(
    Anxiety_Yes = sum(Anxiety == "Yes", na.rm = TRUE),
    Depression_Yes = sum(Depression == "Yes", na.rm = TRUE),
    Behavioral_Problems_Yes = sum(Behavioral_Problems == "Yes", na.rm = TRUE)
  )

# Display table in a simple HTML format using kable
neurodiv_freq_table %>%
  knitr::kable("html", col.names = c("Anxiety", "Depression", "Behavioral Problems"), 
               caption = "Frequency of Mental Health Conditions in Neurodivergent Sample")
Frequency of Mental Health Conditions in Neurodivergent Sample
Anxiety Depression Behavioral Problems
13725 6289 13535

Visualize the distribution of ACE counts among children using a bar plot and calculate cross-tabulations for the number and percentage of participants with Anxiety, Depression, Behavioral Issues, and Flourishing status across ACE categories. Format the table to display values as “n (percentage)” for each ACE category and include a footnote explaining the presentation format.

Expand Code 
# Bar plot of ACE counts categories
ggplot(neurodivergent_data, aes(x = factor(ace_counts), fill = factor(ace_counts))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "ACE Counts", y = "Number of Children", fill = "ACE Counts",
       title = "Distribution of Children by ACE Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))
1
Cross-tabulations (or contingency tables) summarize the relationship between two or more categorical variables. In this case, I’m exploring how the mental health conditions (Anxiety, Depression, Behavioral Issues) and Child Flourishing Index (CFI) are distributed across different levels of ACEs.

Expand Code 
# Define the total number of respondents for each ACE category
ace_totals <- c(
  "0 ACEs" = 12431,  # Total respondents with 0 ACEs
  "1 ACE" = 13433,  # Total respondents with 1 ACE
  "2-3 ACEs" = 12180,  # Total respondents with 2-3 ACEs
  "4+ ACEs" = 6732   # Total respondents with 4+ ACEs
)

# Compute the cross-tabulations, filtering out NA values and excluding "No" responses for Anxiety, Depression, and Behavioral Issues
cross_tabs <- neurodivergent_data %>%
  select(Anxiety, Depression, Behavioral_Problems, CFI_dich, ace_counts) %>%
  pivot_longer(cols = c(Anxiety, Depression, Behavioral_Problems, CFI_dich), names_to = "Condition", values_to = "Status") %>%
  filter(!is.na(Status), !is.na(ace_counts), !(Condition %in% c("Anxiety", "Depression", "Behavioral_Problems") & Status == "No")) %>%  # Filter out NA values and "No" responses
  group_by(Condition, Status, ace_counts) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Condition, Status) %>%
  # Calculate percentage by dividing n by the total number of respondents for each ACE category and create a string in the format n (percent)
  summarise(
    `0 ACEs` = if_else(any(ace_counts == "0 ACEs"), paste0(n[ace_counts == "0 ACEs"], " (", round(n[ace_counts == "0 ACEs"] / ace_totals["0 ACEs"] * 100, 1), ")"), NA_character_),
    `1 ACE` = if_else(any(ace_counts == "1 ACE"), paste0(n[ace_counts == "1 ACE"], " (", round(n[ace_counts == "1 ACE"] / ace_totals["1 ACE"] * 100, 1), ")"), NA_character_),
    `2-3 ACEs` = if_else(any(ace_counts == "2-3 ACEs"), paste0(n[ace_counts == "2-3 ACEs"], " (", round(n[ace_counts == "2-3 ACEs"] / ace_totals["2-3 ACEs"] * 100, 1), ")"), NA_character_),
    `4+ ACEs` = if_else(any(ace_counts == "4+ ACEs"), paste0(n[ace_counts == "4+ ACEs"], " (", round(n[ace_counts == "4+ ACEs"] / ace_totals["4+ ACEs"] * 100, 1), ")"), NA_character_)
  )

# Format the table for presentation using 'kable' with the footnote spanning the full width
cross_tabs %>%
  knitr::kable("html", caption = "Frequency and Percentage of Psychopathology and Flourishing across ACEs Categories") %>%
  kableExtra::kable_styling(full_width = TRUE) %>%  # Ensures the table spans the full width
  kableExtra::footnote(general = "Values are presented as n (percentage) for each ACEs category. Cross-tabulations show the number and percentage of individuals with a specific mental health condition or flourishing status across different ACE categories.",
                       footnote_as_chunk = TRUE)  # Formats the footnote to span full width
Frequency and Percentage of Psychopathology and Flourishing across ACEs Categories
Condition Status 0 ACEs 1 ACE 2-3 ACEs 4+ ACEs
Anxiety Yes 2695 (21.7) 3549 (26.4) 4257 (35) 3224 (47.9)
Behavioral_Problems Yes 2530 (20.4) 3561 (26.5) 4189 (34.4) 3255 (48.4)
CFI_dich Not Flourishing 11140 (89.6) 12358 (92) 11492 (94.4) 6521 (96.9)
CFI_dich Flourishing 1291 (10.4) 1075 (8) 688 (5.6) 211 (3.1)
Depression Yes 804 (6.5) 1222 (9.1) 2112 (17.3) 2151 (32)
Note: Values are presented as n (percentage) for each ACEs category. Cross-tabulations show the number and percentage of individuals with a specific mental health condition or flourishing status across different ACE categories.

Inferential Statistics

Run three separate logistic regression models to examine the effects of ACE counts, race, sex, socioeconomic status, caregiver education, and child age on Anxiety, Depression, and Behavioral Issues outcomes. Present the exponentiated coefficients (odds ratios) and bold significant p-values (p < 0.05). Merge the summary tables of all three models into one, with each model labeled appropriately for clear comparison.

Expand Code 
# Logistic regression models without interactions
model_anx <- glm(Anxiety ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

model_dep <- glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

model_beh <- glm(Behavioral_Problems ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years,
               data = neurodivergent_data, family = binomial)

# Create a summary table for the Anixety model
tbl_anx <- tbl_regression(model_anx, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Depression model
tbl_dep <- tbl_regression(model_dep, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Behavioral Problems model
tbl_beh <- tbl_regression(model_beh, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Merge all tables into one
tbl_merged <- tbl_merge(
  tbls = list(tbl_anx, tbl_dep, tbl_beh),
  tab_spanner = c("Anixety Model", "Depression Model", "Behavioral Problems Model")
)

# Display the combined table
tbl_merged
Characteristic
Anixety Model
Depression Model
Behavioral Problems Model
OR1 95% CI1 p-value OR1 95% CI1 p-value OR1 95% CI1 p-value
ACEs








    0 ACEs — —
— —
— —
    1 ACE 1.30 1.10, 1.53 0.002 1.23 1.01, 1.49 0.042 1.14 1.00, 1.31 0.054
    2-3 ACEs 1.49 1.26, 1.76 <0.001 1.45 1.20, 1.74 <0.001 1.11 0.97, 1.27 0.13
    4+ ACEs 1.86 1.53, 2.26 <0.001 1.78 1.46, 2.18 <0.001 1.37 1.18, 1.59 <0.001
Child Race








    Non-Hispanic White — —
— —
— —
    Non-Hispanic Black or African American 0.77 0.59, 1.02 0.057 0.94 0.72, 1.25 0.7 1.03 0.86, 1.24 0.7
    Non-Hispanic American Indian or Alaska Native 0.52 0.31, 0.93 0.019 1.35 0.68, 3.06 0.4 0.91 0.57, 1.52 0.7
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander 0.63 0.43, 0.95 0.021 1.34 0.83, 2.30 0.3 1.14 0.84, 1.58 0.4
    Non-Hispanic Other or Mixed Race 0.89 0.71, 1.14 0.4 1.01 0.80, 1.29 >0.9 0.93 0.78, 1.10 0.4
    Hispanic/Latino of Any Race 0.71 0.59, 0.85 <0.001 0.84 0.69, 1.02 0.070 0.90 0.78, 1.05 0.2
Child Sex








    Male — —
— —
— —
    Female 1.34 1.18, 1.51 <0.001 1.42 1.25, 1.62 <0.001 1.21 1.09, 1.35 <0.001
Federal Poverty Level








    <100% — —
— —
— —
    100%-199% 0.94 0.75, 1.17 0.6 0.89 0.71, 1.12 0.3 0.91 0.77, 1.08 0.3
    200%-299% 1.03 0.82, 1.31 0.8 0.84 0.66, 1.06 0.14 0.83 0.70, 0.99 0.043
    300%-399% 0.82 0.64, 1.04 0.10 0.79 0.61, 1.01 0.062 0.68 0.57, 0.82 <0.001
    ≥400% 0.91 0.72, 1.14 0.4 0.69 0.54, 0.87 0.002 0.71 0.59, 0.84 <0.001
Caregiver Education








    Less Than High School — —
— —
— —
    High School/Vocational 1.18 0.77, 1.77 0.4 1.10 0.71, 1.64 0.7 1.01 0.73, 1.37 >0.9
    Some College/Associate Degree 1.18 0.77, 1.74 0.4 1.13 0.74, 1.68 0.6 0.99 0.72, 1.34 >0.9
    College Degree/Higher 1.30 0.85, 1.93 0.2 1.03 0.67, 1.53 >0.9 0.87 0.63, 1.17 0.4
Child Age 0.97 0.95, 0.99 <0.001 0.96 0.93, 0.98 <0.001 0.86 0.84, 0.87 <0.001
1 OR = Odds Ratio, CI = Confidence Interval

Generate a combined plot showing the odds ratios from three separate logistic regression models (Anxiety, Depression, and Behavioral Issues). Configure the figure dimensions to 8x8, with confidence intervals (95%) and p-values displayed. Use different markers and colors for each model, set axis limits between 0.2 and 5, and include a legend labeling each model accordingly.

Expand Code 
# Combine the plots into a list vector
all.models <- list()
all.models[[1]] <- model_anx
all.models[[2]] <- model_dep
all.models[[3]] <- model_beh

# Plot the combined version
plot_models(all.models, axis.lim = c(0.2, 5), dot.size = 1, show.p = TRUE, ci.lvl = 0.95, title = "Odds Ratios by Model", legend.title = "Model", m.labels = c("Anxiety", "Depression", "Behavioral Problems"))

Run logistic regression models that incorporate interaction terms between the ACE counts and the Child Flourishing Index (CFI_dich) for three different mental health outcomes (Anxiety, Depression, Behavioral Issues). Then, create individual summary tables for each model (without exponentiating the coefficients), highlighting statistically significant p-values. Finally, merge the tables into a single output, displaying all models side by side with appropriate labels and formatting.

Expand Code 
# Logistic regression model with CFI_dich as a moderator for Anxiety
model_anx_mod <- glm(Anxiety ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)

# Create a summary table for the Anixety model
tbl_anx_mod <- tbl_regression(model_anx_mod, exponentiate = FALSE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Logistic regression model with CFI_dich as a moderator for Depression
model_dep_mod <- glm(Depression ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)

# Create a summary table for the Depression model
tbl_dep_mod <- tbl_regression(model_dep_mod, exponentiate = FALSE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Logistic regression model with CFI_dich as a moderator for Behavioral Problems
model_beh_mod <- glm(Behavioral_Problems ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = neurodivergent_data, family = binomial)

# Create a summary table for the Behavioral Problems model
tbl_beh_mod <- tbl_regression(model_beh_mod, exponentiate = FALSE, label = list( 
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Merge all tables into one
mod_tbl_merged <- tbl_merge(
  tbls = list(tbl_anx_mod, tbl_dep_mod, tbl_beh_mod),
  tab_spanner = c("Anixety Model", "Depression Model", "Behavioral Problems Model")
)

# Display the combined table
mod_tbl_merged
1
The interaction anxiety model did not yield significant findings, so we do not proceed with post-hoc analyses.
2
The interaction depression model did yield significant findinds, so we do proceed with post-hoc analyses.
3
The interaction behavioral problems model did yield significant findinds, so we do proceed with post-hoc analyses.
Characteristic
Anixety Model
Depression Model
Behavioral Problems Model
log(OR)1 95% CI1 p-value log(OR)1 95% CI1 p-value log(OR)1 95% CI1 p-value
ACEs








    0 ACEs — —
— —
— —
    1 ACE 0.26 0.09, 0.43 0.003 0.19 -0.01, 0.39 0.065 0.14 0.00, 0.28 0.043
    2-3 ACEs 0.38 0.20, 0.55 <0.001 0.32 0.13, 0.51 <0.001 0.11 -0.03, 0.25 0.12
    4+ ACEs 0.58 0.38, 0.78 <0.001 0.54 0.33, 0.74 <0.001 0.31 0.16, 0.46 <0.001
Child Flourishing Index








    Not Flourishing — —
— —
— —
    Flourishing -0.99 -1.4, -0.50 <0.001 -1.4 -2.2, -0.61 <0.001 -0.71 -1.3, -0.13 0.013
Child Race








    Non-Hispanic White — —
— —
— —
    Non-Hispanic Black or African American -0.27 -0.54, 0.01 0.050 -0.06 -0.33, 0.23 0.7 0.03 -0.15, 0.22 0.7
    Non-Hispanic American Indian or Alaska Native -0.66 -1.2, -0.07 0.018 0.28 -0.40, 1.1 0.5 -0.11 -0.57, 0.41 0.7
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander -0.48 -0.85, -0.07 0.016 0.29 -0.19, 0.83 0.3 0.11 -0.20, 0.44 0.5
    Non-Hispanic Other or Mixed Race -0.13 -0.36, 0.12 0.3 0.01 -0.23, 0.25 >0.9 -0.07 -0.24, 0.10 0.4
    Hispanic/Latino of Any Race -0.34 -0.52, -0.16 <0.001 -0.18 -0.37, 0.02 0.075 -0.10 -0.25, 0.05 0.2
Child Sex








    Male — —
— —
— —
    Female 0.30 0.18, 0.42 <0.001 0.36 0.23, 0.49 <0.001 0.19 0.09, 0.30 <0.001
Federal Poverty Level








    <100% — —
— —
— —
    100%-199% -0.06 -0.28, 0.16 0.6 -0.11 -0.34, 0.12 0.4 -0.11 -0.29, 0.06 0.2
    200%-299% 0.03 -0.20, 0.27 0.8 -0.18 -0.42, 0.06 0.14 -0.20 -0.38, -0.02 0.026
    300%-399% -0.19 -0.43, 0.05 0.12 -0.24 -0.50, 0.01 0.062 -0.40 -0.59, -0.21 <0.001
    ≥400% -0.09 -0.32, 0.14 0.5 -0.37 -0.61, -0.14 0.002 -0.37 -0.55, -0.19 <0.001
Caregiver Education








    Less Than High School — —
— —
— —
    High School/Vocational 0.14 -0.29, 0.55 0.5 0.09 -0.34, 0.50 0.7 -0.01 -0.34, 0.30 >0.9
    Some College/Associate Degree 0.13 -0.30, 0.52 0.5 0.12 -0.30, 0.51 0.6 -0.04 -0.36, 0.27 0.8
    College Degree/Higher 0.22 -0.21, 0.61 0.3 0.02 -0.41, 0.42 >0.9 -0.17 -0.49, 0.14 0.3
Child Age -0.03 -0.05, -0.01 0.002 -0.04 -0.07, -0.02 <0.001 -0.15 -0.17, -0.14 <0.001
ACEs * Child Flourishing Index








    1 ACE * Flourishing -0.23 -0.91, 0.46 0.5 -0.13 -1.3, 1.0 0.8 -0.52 -1.3, 0.24 0.2
    2-3 ACEs * Flourishing 0.10 -0.58, 0.80 0.8 1.3 0.27, 2.5 0.017 -0.54 -1.3, 0.20 0.2
    4+ ACEs * Flourishing 0.54 -0.44, 1.7 0.3 1.0 -0.08, 2.2 0.077 -1.1 -2.1, -0.19 0.019
1 OR = Odds Ratio, CI = Confidence Interval

Subset the data into flourishing and non-flourishing groups. Then, run logistic regression models for Depression and Behavioral Problems separately for each group, adjusting for covariates (ACE counts, race, sex, SES, education, and age). Create summary tables for each model, with p-values highlighted for significance. Finally, merge and display the tables for both flourishing and non-flourishing samples, stratified by Depression and Behavioral Problems.

Expand Code 
# Subset the flourishing sample
flourishing_sample <- neurodivergent_data %>%
  filter(CFI_dich == "Flourishing")

# Subset the not flourishing sample
not_flourishing_sample <- neurodivergent_data %>%
  filter(CFI_dich == "Not Flourishing")

# Stratified logistic regression models
flourish_model_dep <- glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = flourishing_sample, family = binomial)

non_flourish_model_dep <- glm(Depression ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = not_flourishing_sample, family = binomial)

flourish_model_beh <- glm(Behavioral_Problems ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = flourishing_sample, family = binomial)

non_flourish_model_beh <- glm(Behavioral_Problems ~ ace_counts + sc_race_cat + sc_sex_cat + fpl_cat + highgrade_tvis_cat + sc_age_years, data = not_flourishing_sample, family = binomial)

# Create a summary table for the Flourishing Depression model
tbl_dep_flourish <- tbl_regression(flourish_model_dep, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Non-Flourishing Depression model
tbl_dep_non_flourish <- tbl_regression(non_flourish_model_dep, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Flourishing Behavioral Problems model
tbl_beh_flourish <- tbl_regression(flourish_model_beh, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Create a summary table for the Non-Flourishing Behavioral Problems model
tbl_beh_non_flourish <- tbl_regression(non_flourish_model_beh, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      fpl_cat ~ "Federal Poverty Level",
      highgrade_tvis_cat ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex"
    )) |>
  bold_p(t = 0.05)

# Merge the flourishing model tables into one
flourish_tbl_combined <- tbl_merge(
  tbls = list(tbl_dep_flourish, tbl_beh_flourish),
  tab_spanner = c("Depression Model", "Behavioral Problems Model")
)

# Display the combined flourishing model tables
flourish_tbl_combined
Characteristic
Depression Model
Behavioral Problems Model
OR1 95% CI1 p-value OR1 95% CI1 p-value
ACEs





    0 ACEs — —
— —
    1 ACE 1.18 0.34, 4.11 0.8 0.71 0.32, 1.57 0.4
    2-3 ACEs 7.81 2.37, 28.8 0.001 0.64 0.29, 1.38 0.3
    4+ ACEs 8.99 2.36, 38.1 0.002 0.36 0.12, 1.05 0.064
Child Race





    Non-Hispanic White — —
— —
    Non-Hispanic Black or African American 0.42 0.09, 2.09 0.3 0.80 0.29, 2.28 0.7
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander 20,679,597 0.00,
 >0.9 3,087,910 0.00,
 >0.9
    Non-Hispanic Other or Mixed Race 0.68 0.10, 6.07 0.7 1.42 0.47, 4.94 0.6
    Hispanic/Latino of Any Race 1.03 0.29, 4.22 >0.9 0.40 0.16, 0.96 0.042
    Non-Hispanic American Indian or Alaska Native


0.20 0.01, 5.72 0.3
Child Sex





    Male — —
— —
    Female 1.12 0.48, 2.63 0.8 1.72 0.88, 3.48 0.12
Federal Poverty Level





    <100% — —
— —
    100%-199% 1.65 0.41, 7.05 0.5 0.73 0.27, 1.93 0.5
    200%-299% 1.60 0.33, 8.40 0.6 0.99 0.35, 2.75 >0.9
    300%-399% 0.84 0.18, 4.05 0.8 0.93 0.30, 2.87 >0.9
    ≥400% 2.62 0.58, 12.5 0.2 0.68 0.23, 1.92 0.5
Caregiver Education





    Less Than High School — —
— —
    High School/Vocational 0.94 0.04, 8.79 >0.9 0.90 0.20, 3.61 0.9
    Some College/Associate Degree 0.44 0.02, 3.64 0.5 1.00 0.22, 3.95 >0.9
    College Degree/Higher 0.33 0.01, 2.97 0.4 0.62 0.13, 2.53 0.5
Child Age 0.91 0.75, 1.10 0.3 0.86 0.78, 0.95 0.003
1 OR = Odds Ratio, CI = Confidence Interval
Expand Code 
# Merge the non-flourishing model tables into one
non_flourish_tbl_combined <- tbl_merge(
  tbls = list(tbl_dep_non_flourish, tbl_beh_non_flourish),
  tab_spanner = c("Depression Model", "Behavioral Problems Model")
)

# Display the combined non-flourishing model tables
non_flourish_tbl_combined
Characteristic
Depression Model
Behavioral Problems Model
OR1 95% CI1 p-value OR1 95% CI1 p-value
ACEs





    0 ACEs — —
— —
    1 ACE 1.21 0.99, 1.48 0.067 1.15 1.00, 1.33 0.042
    2-3 ACEs 1.38 1.14, 1.66 0.001 1.12 0.97, 1.28 0.11
    4+ ACEs 1.70 1.38, 2.09 <0.001 1.37 1.17, 1.59 <0.001
Child Race





    Non-Hispanic White — —
— —
    Non-Hispanic Black or African American 0.97 0.73, 1.29 0.8 1.04 0.87, 1.26 0.7
    Non-Hispanic American Indian or Alaska Native 1.33 0.67, 3.02 0.5 0.93 0.58, 1.58 0.8
    Non-Hispanic Asian, Native Hawaiian, or Pacific Islander 1.30 0.80, 2.24 0.3 1.10 0.81, 1.53 0.6
    Non-Hispanic Other or Mixed Race 1.01 0.80, 1.30 >0.9 0.92 0.77, 1.10 0.3
    Hispanic/Latino of Any Race 0.83 0.69, 1.02 0.067 0.92 0.80, 1.07 0.3
Child Sex





    Male — —
— —
    Female 1.44 1.26, 1.63 <0.001 1.20 1.08, 1.34 <0.001
Federal Poverty Level





    <100% — —
— —
    100%-199% 0.88 0.70, 1.11 0.3 0.89 0.75, 1.06 0.2
    200%-299% 0.83 0.65, 1.05 0.13 0.81 0.68, 0.97 0.024
    300%-399% 0.78 0.60, 1.01 0.059 0.66 0.55, 0.80 <0.001
    ≥400% 0.67 0.53, 0.85 <0.001 0.69 0.58, 0.83 <0.001
Caregiver Education





    Less Than High School — —
— —
    High School/Vocational 1.11 0.72, 1.67 0.6 1.00 0.71, 1.37 >0.9
    Some College/Associate Degree 1.16 0.75, 1.73 0.5 0.97 0.70, 1.33 0.9
    College Degree/Higher 1.05 0.68, 1.57 0.8 0.85 0.61, 1.17 0.3
Child Age 0.96 0.94, 0.98 <0.001 0.86 0.84, 0.87 <0.001
1 OR = Odds Ratio, CI = Confidence Interval

Generate odds ratio plots for Depression and Behavioral Problems stratified by flourishing status. First, create separate plots for the flourishing group and the non-flourishing group, adjusting the x-axis limits, confidence intervals, and dot size. The plots are combined for each group to visualize the odds ratios for the Depression and Behavioral Problems models, and labeled accordingly for each condition within both flourishing and non-flourishing samples.

Expand Code 
# Combine the flourishing plots into a list vector
all.flourish.models <- list()
all.flourish.models[[1]] <- flourish_model_dep
all.flourish.models[[2]] <- flourish_model_beh

# Plot the combined version
plot_models(all.flourish.models, axis.lim = c(0.2, 5), dot.size = 1, show.p = TRUE, ci.lvl = 0.95, title = "Odds Ratios by Flourishing Stratified Model", legend.title = "Model", m.labels = c("Depression", "Behavioral Problems"))

Expand Code 
# Combine the non-flourishing plots into a list vector
all.non.flourish.models <- list()
all.non.flourish.models[[1]] <- non_flourish_model_dep
all.non.flourish.models[[2]] <- non_flourish_model_beh

# Plot the combined version
plot_models(all.flourish.models, axis.lim = c(0.2, 5), dot.size = 1, show.p = TRUE, ci.lvl = 0.95, title = "Odds Ratios by Non-Flourishing Stratified Model", legend.title = "Model", m.labels = c("Depression", "Behavioral Problems"))

Expand Code 
# Create a diagram representing the analysis flow
diagram <- grViz("
  digraph flowchart {
  
    # Node definitions
    node [shape = box, style = filled, fillcolor = lightgrey, fontname = Helvetica]
    X [label = 'Adverse Childhood Experiences', width = 2]
    Y [label = 'Anxiety, Depression, & Behavioral Problems', width = 2]
    W [label = 'Child Flourishing', shape = box, style = filled, fillcolor = lightgrey, fontname = Helvetica]
    Interaction [label = 'ACEs * Child Flourishing Interaction', shape = box, style = filled, fillcolor = lightgrey, fontname = Helvetica]

    # Define edge connections
    X -> Y [label = 'a: ##']
    W -> Y [label = 'b: ##']
    Interaction -> Y [label = 'c: ##']  # Interaction effect
    X -> Interaction [label = '', style = invis]
    W -> Interaction [label = '', style = invis]
  
    # Graph layout
    rankdir = LR;  # Left to right layout
  }
")

# Render the diagram
diagram

Session Information

To enhance reproducibility, details about the working environment used for these analyses can be found below.

Expand Code 
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS 15.0.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DiagrammeR_1.0.11    vtable_1.4.6         sjmisc_2.8.10       
 [4] sjlabelled_1.2.0     patchwork_1.2.0      broom.helpers_1.17.0
 [7] ggstats_0.7.0        gtsummary_2.0.3      margins_0.3.28      
[10] sjPlot_2.8.16        webshot2_0.1.1       kableExtra_1.4.0    
[13] broom_1.0.6          Hmisc_5.1-2          knitr_1.48          
[16] lubridate_1.9.3      forcats_1.0.0        stringr_1.5.1       
[19] purrr_1.0.2          readr_2.1.5          tibble_3.2.1        
[22] ggplot2_3.5.1        tidyverse_2.0.0      tidyr_1.3.1         
[25] dplyr_1.1.4          pacman_0.5.1        

loaded via a namespace (and not attached):
 [1] gridExtra_2.3       rlang_1.1.4         magrittr_2.0.3     
 [4] snakecase_0.11.1    prediction_0.3.18   compiler_4.4.1     
 [7] systemfonts_1.1.0   vctrs_0.6.5         pkgconfig_2.0.3    
[10] crayon_1.5.3        fastmap_1.2.0       backports_1.5.0    
[13] labeling_0.4.3      effectsize_0.8.9    utf8_1.2.4         
[16] PupillometryR_0.0.5 promises_1.3.0      rmarkdown_2.27     
[19] markdown_1.13       tzdb_0.4.0          haven_2.5.4        
[22] ps_1.7.6            bit_4.0.5           xfun_0.45          
[25] labelled_2.13.0     jsonlite_1.8.8      highr_0.11         
[28] later_1.3.2         ggeffects_1.7.0     parallel_4.4.1     
[31] cluster_2.1.6       R6_2.5.1            stringi_1.8.4      
[34] RColorBrewer_1.1-3  rpart_4.1.23        Rcpp_1.0.12        
[37] parameters_0.22.0   base64enc_0.1-3     nnet_7.3-19        
[40] timechange_0.3.0    tidyselect_1.2.1    rstudioapi_0.16.0  
[43] yaml_2.3.9          websocket_1.4.1     processx_3.8.4     
[46] bayestestR_0.13.2   withr_3.0.0         evaluate_0.24.0    
[49] foreign_0.8-86      xml2_1.3.6          pillar_1.9.0       
[52] checkmate_2.3.1     insight_0.20.1      generics_0.1.3     
[55] vroom_1.6.5         chromote_0.2.0      hms_1.1.3          
[58] commonmark_1.9.1    munsell_0.5.1       scales_1.3.0       
[61] glue_1.8.0          tools_4.4.1         data.table_1.15.4  
[64] visNetwork_2.1.2    grid_4.4.1          cards_0.3.0        
[67] datawizard_0.11.0   colorspace_2.1-0    performance_0.12.0 
[70] htmlTable_2.4.2     Formula_1.2-5       cli_3.6.3          
[73] fansi_1.0.6         viridisLite_0.4.2   gt_0.11.1          
[76] svglite_2.1.3       sjstats_0.19.0      gtable_0.3.5       
[79] sass_0.4.9          digest_0.6.36       htmlwidgets_1.6.4  
[82] farver_2.1.2        htmltools_0.5.8.1   lifecycle_1.0.4    
[85] bit64_4.0.5         MASS_7.3-60.2

This document was made by Adrián Medina as part of the Slopen Laboratory, based on the Lumo quarto template.